%% 清空环境
clear; clc; close all;

%% 固定参数
a0 = -2;
a2 = 0.14;
b0 = 1.3;
b1 = 0.1;
b2 = 1;
d1 = 0;
d2 = 0;

% 设定 a1 扫描范围
a1_min = 0;
a1_max = 6*pi;
num_a1 = 601;  % a1 取样点数
a1_values = linspace(a1_min, a1_max, num_a1);

% 迭代设置
N_total = 1000;     % 总迭代步数
N_trans = 200;     % 舍弃前 200 步
sample_interval = 10;  % 分叉图采样间隔

%% 预分配存储
LE_result = zeros(num_a1, 3);   % 用于保存 (LE1, LE2, LE3)
a1_bif = [];     % 分叉图 横坐标 a1
x_bif  = [];     % 分叉图 纵坐标 x

%% 扫描 a1
for i = 1:num_a1
    a1 = a1_values(i);
    
    % 初始状态
    x = 0.1;
    y = 0.1;
    z = 0.1;
    
    % 李指数相关
    Q = eye(3);          % QR法中的正交矩阵
    sum_log = zeros(1,3);% 累计 log|Rii|
    
    % 迭代计算
    count_bif = 0;  % 用于记录分叉数据
    
    for n = 1:N_total
        % 映射更新
        x_new = x + y;
        y_new = sin(a0 * y^2 * sin(a1 + a2*x)) + b2*sin(z) + d1;
        z_new = b0*y + b1*z + d2;
        
        % 计算 Jacobian 矩阵
        % ----------------------------------
        %   J = [ 1,     1,          0
        %         M1,    M2,   b2*cos(z)
        %         0,     b0,        b1  ]
        % 其中:
        %   M1 = a0*a2*y^2*cos(a1+a2*x)*cos(a0*y^2*sin(a1+a2*x))
        %   M2 = 2*a0*y*cos(a0*y^2*sin(a1+a2*x))*sin(a1+a2*x)
        
        theta = a1 + a2*x;
        sin_theta = sin(theta);
        cos_theta = cos(theta);
        arg = a0 * y^2 * sin_theta;
        cos_arg = cos(arg);
        
        M1 = a0 * a2 * y^2 * cos_theta * cos_arg;
        M2 = 2 * a0 * y * cos_arg * sin_theta;
        
        J = [1, 1, 0;
             M1, M2, b2*cos(z);
             0, b0, b1];
        
        % 更新 Q 矩阵 (QR 正交化)
        Q = J * Q;
        [Q, R] = qr(Q);
        
        % 累计 log|Rii|
        if n > N_trans
            diagR = diag(R);
            sum_log = sum_log + log(abs(diagR))';
        end
        
        % 更新状态
        x = x_new;
        y = y_new;
        z = z_new;
        
        % 收集分叉图数据 (舍弃瞬态后每隔 sample_interval 记录一次)
        if (n > N_trans) && (mod(n, sample_interval) == 0)
            x_bif(end+1) = x;
            a1_bif(end+1) = a1;
            count_bif = count_bif + 1;
        end
    end
    
    % 计算该 a1 下的 Lyapunov 指数
    num_effective = N_total - N_trans;
    LE = sum_log / num_effective;
    LE = sort(LE, 'descend');  % 降序排列 (LE1 >= LE2 >= LE3)
    LE_result(i,:) = LE;
    
    fprintf('Progress: %.1f%% \n', (i/num_a1)*100 );
end

%% 绘图
figure;

% (1) 上子图：Lyapunov 指数
subplot(2,1,1);
hold on; grid on;
plot(a1_values, LE_result(:,1), 'r-', 'LineWidth',1, 'DisplayName','LE_1');
plot(a1_values, LE_result(:,2), 'b-', 'LineWidth',1, 'DisplayName','LE_2');
plot(a1_values, LE_result(:,3), 'g-', 'LineWidth',1, 'DisplayName','LE_3');
plot(a1_values, zeros(size(a1_values)), 'k--', 'LineWidth',1);
hold off;
xlabel('a_1','FontSize',12);
ylabel('LE','FontSize',12);
title('a_1 李指数图','FontSize',14);
legend('Location','best');

% (2) 下子图：分叉图 (x_n vs. a_1)
subplot(2,1,2);
plot(a1_bif, x_bif, 'r.', 'MarkerSize',1);
grid on;
xlabel('a_1','FontSize',12);
ylabel('x_n','FontSize',12);
title(' x -. a_1 分岔图','FontSize',14);
